#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

#define PRECISION 3500  // 设置精度（以比特为单位），至少需要 3320 位来表示 1000 位十进制数字

void calculate_pi(mpf_t pi) {
    // 初始化变量，设置所需的精度
    mpf_t a, b, t, p, a_next, b_next, temp1, temp2, temp3;
    mpf_init2(a, PRECISION);
    mpf_init2(b, PRECISION);
    mpf_init2(t, PRECISION);
    mpf_init2(p, PRECISION);
    mpf_init2(a_next, PRECISION);
    mpf_init2(b_next, PRECISION);
    mpf_init2(temp1, PRECISION);
    mpf_init2(temp2, PRECISION);
    mpf_init2(temp3, PRECISION);

    // 设置初始值
    mpf_set_d(a, 1.0);

    // b = 1 / sqrt(2)
    mpf_sqrt_ui(b, 2);       // b = sqrt(2)
    mpf_ui_div(b, 1, b);     // b = 1 / sqrt(2)

    mpf_set_d(t, 0.25);      // t = 1/4
    mpf_set_d(p, 1.0);       // p = 1

    int N = 11;  // 迭代次数，11 次迭代可以得到超过 1000 位精度

    for (int i = 0; i < N; i++) {
        // a_next = (a + b) / 2
        mpf_add(a_next, a, b);
        mpf_div_ui(a_next, a_next, 2);

        // b_next = sqrt(a * b)
        mpf_mul(b_next, a, b);
        mpf_sqrt(b_next, b_next);

        // temp1 = a - a_next
        mpf_sub(temp1, a, a_next);

        // temp2 = (a - a_next)^2
        mpf_mul(temp2, temp1, temp1);

        // temp3 = p * temp2
        mpf_mul(temp3, p, temp2);

        // t = t - temp3
        mpf_sub(t, t, temp3);

        // 更新 p：p = 2 * p
        mpf_mul_ui(p, p, 2);

        // 更新 a 和 b
        mpf_set(a, a_next);
        mpf_set(b, b_next);
    }

    // 计算 PI：pi = ((a + b)^2) / (4 * t)
    mpf_add(temp1, a, b);        // temp1 = a + b
    mpf_mul(temp1, temp1, temp1); // temp1 = (a + b)^2
    mpf_mul_ui(temp2, t, 4);      // temp2 = 4 * t
    mpf_div(pi, temp1, temp2);    // pi = temp1 / temp2

    // 清理内存
    mpf_clear(a);
    mpf_clear(b);
    mpf_clear(t);
    mpf_clear(p);
    mpf_clear(a_next);
    mpf_clear(b_next);
    mpf_clear(temp1);
    mpf_clear(temp2);
    mpf_clear(temp3);
}

int main() {
    mpf_t pi;
    mpf_init2(pi, PRECISION); // 设置精度

    calculate_pi(pi);

    // 打印 PI，保留小数点后 1000 位
    gmp_printf("PI = %.1000Ff\n", pi);

    mpf_clear(pi); // 清理内存
    return 0;
}

